suppressPackageStartupMessages({
library("here")
library("tidyverse")
library("MplusAutomation")
library("gt")
library("glue")
library("kableExtra")
library("misty")
library("lavaan")
library("AICcmodavg")
library("nonnest2")
library("DiagrammeR")
library("lavaan")
library("tidyLPA")
library("semTools")
library("brms")
library("MBESS")
library("ufs")
library("robmed")
library("careless")
library("mice")
library("psych")
library("BayesFactor")
library("effectsize")
library("tidybayes")
library("emmeans")
library("bayesplot")
library("patchwork")
library("bmlm")
})
# install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
options("max.print" = .Machine$integer.max)
# Make random things reproducible
set.seed(1234)
options(mc.cores = 4)
bayesplot_theme_set()
source(here::here("src", "functions", "funs_add_neoffi60_subscales.R"))
source(here::here("src", "functions", "funs_correct_iesr_scores.R"))
source(here::here("src", "functions", "funs_plot_job_qualification.R"))
source(here::here("src", "functions", "funs_generate_all_items_df.R"))
scale_this <- function(x) as.vector(scale(x))
sum_coding <- function(x, lvls = levels(x)) {
# codes the first category with -1
nlvls <- length(lvls)
stopifnot(nlvls > 1)
cont <- diag(nlvls)[, -nlvls, drop = FALSE]
cont[nlvls, ] <- -1
cont <- cont[c(nlvls, 1:(nlvls - 1)), , drop = FALSE]
colnames(cont) <- lvls[-1]
x <- factor(x, levels = lvls)
contrasts(x) <- cont
x
}
all_items <- generate_all_items_df()
There is a problem with IES-R, in the control group. I shift the control distribution of IES-R towards lower values.
```r
temp <- correct_iesr_scores(all_items)
all_items <- temp
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuZ2dwbG90KGFsbF9pdGVtcywgYWVzKHggPSBpZXNyX3RzLCBjb2xvdXIgPSBpc19yZXNjdWVfd29ya2VyKSkgK1xuICBnZW9tX2RlbnNpdHkoKVxuYGBgIn0= -->
```r
ggplot(all_items, aes(x = iesr_ts, colour = is_rescue_worker)) +
geom_density()
all_items |>
ggplot(aes(y=iesr_ts, fill=is_rescue_worker)) +
geom_boxplot() +
scale_y_log10()
all_items$y <- all_items$iesr_ts + 0.001
all_items |>
group_by(is_rescue_worker) |>
summarize(
# gm_iesr = geometric.mean(y, na.rm = TRUE),
avg_iesr = mean(y, na.rm = TRUE)
)
There are some cases in which the RWs of Toscana and Lombardia did not have the proper qualification. They are assigned to the category of RWs.
all_items$idx <- 1:nrow(all_items)
all_items$is_rescue_worker <- as.character(all_items$is_rescue_worker)
all_items$is_rescue_worker <- ifelse(
all_items$idx < 746, "Si", all_items$is_rescue_worker
)
all_items$idx <- NULL
all_items$commeetee_location <- ifelse(
all_items$is_rescue_worker == "No", "None", all_items$red_cross_commeetee_location
) |>
factor()
all_items$commeetee_location <- ifelse(
all_items$is_rescue_worker == "No" & all_items$age < 25, "students", all_items$commeetee_location
)
all_items$commeetee <- ifelse(
all_items$is_rescue_worker == "No" & all_items$age >= 25, "community_sample",
all_items$commeetee_location
)
all_items$commeetee <-
ifelse(is.na(all_items$commeetee), "community_sample", all_items$commeetee) |>
factor()
sum(is.na(all_items$commeetee))
[1] 0
```r
```r
```r
# all_items <- all_items %>%
# mutate(job_qualification = case_when(
# is_rescue_worker==\Si\ & job_qualification == \non_rescue_worker\ ~ \team_member\,
# TRUE ~ job_qualification))
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
## IES-R as a function of group
iesr_ts | trunc(lb = 0) ~ is_rescue_worker + (1 | commeetee),
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubTAgPC0gYnJtKFxuICBiZihpZXNyX3RzIH4gaXNfcmVzY3VlX3dvcmtlciksXG4gIGZhbWlseSA9IGh1cmRsZV9nYW1tYSgpLFxuICBkYXRhID0gYWxsX2l0ZW1zLFxuICBiYWNrZW5kID0gXCJjbWRzdGFuclwiXG4gICMgYWxnb3JpdGhtID0gXCJtZWFuZmllbGRcIlxuKVxuYGBgIn0= -->
```r
m0 <- brm(
bf(iesr_ts ~ is_rescue_worker),
family = hurdle_gamma(),
data = all_items,
backend = "cmdstanr"
# algorithm = "meanfield"
)
pp_check(m0)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m0)
Family: hurdle_gamma
Links: mu = log; shape = identity; hu = identity
Formula: iesr_ts ~ is_rescue_worker
Data: all_items (Number of observations: 1068)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.69 0.05 2.60 2.79 1.00 4431 2990
is_rescue_workerSi 0.26 0.06 0.15 0.37 1.00 4401 3236
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 1.41 0.06 1.30 1.52 1.00 4553 2820
hu 0.05 0.01 0.04 0.06 1.00 4772 3170
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m0, "is_rescue_worker"
)
plot(me, points = FALSE)
BFt <- BayesFactor::ttestBF(
all_items$ies_ts[all_items$is_rescue_worker == "Si"],
all_items$ies_ts[all_items$is_rescue_worker == "No"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
------------------------
0.28 | [0.15, 0.41]
- Estimated using pooled SD.
Supported families are: ‘acat’, ‘asym_laplace’, ‘bernoulli’, ‘beta’, ‘beta_binomial’, ‘binomial’, ‘categorical’, ‘com_poisson’, ‘cox’, ‘cratio’, ‘cumulative’, ‘custom’, ‘dirichlet’, ‘dirichlet2’, ‘discrete_weibull’, ‘exgaussian’, ‘exponential’, ‘frechet’, ‘gamma’, ‘gaussian’, ‘gen_extreme_value’, ‘geometric’, ‘hurdle_cumulative’, ‘hurdle_gamma’, ‘hurdle_lognormal’, ‘hurdle_negbinomial’, ‘hurdle_poisson’, ‘info’, ‘inverse.gaussian’, ‘logistic_normal’, ‘lognormal’, ‘multinomial’, ‘negbinomial’, ‘negbinomial2’, ‘poisson’, ‘shifted_lognormal’, ‘skew_normal’, ‘sratio’, ‘student’, ‘von_mises’, ‘weibull’, ‘wiener’, ‘zero_inflated_asym_laplace’, ‘zero_inflated_beta’, ‘zero_inflated_beta_binomial’, ‘zero_inflated_binomial’, ‘zero_inflated_negbinomial’, ‘zero_inflated_poisson’, ‘zero_one_inflated_beta’
The sk, ch, mi sub-scales are coded so that high values indicate high self-compassion levels. The sj, is, oi sub-scales are coded so that high values indicate low self-compassion levels.
The ts_sc score has been computed by reversing the coding of the items of the sj, is, oi sub-scales (so that they indicate the absence of self-judgment, absence of isolation, absence of over-identification).
scs_subscales <- with(all_items, data.frame(sk, ch, mi, sj, is, oi, scs_ts))
cor(scs_subscales) |> round(2)
sk ch mi sj is oi scs_ts
sk 1.00 0.52 0.58 -0.39 -0.28 -0.24 0.71
ch 0.52 1.00 0.49 -0.01 -0.03 -0.04 0.45
mi 0.58 0.49 1.00 -0.19 -0.33 -0.35 0.66
sj -0.39 -0.01 -0.19 1.00 0.67 0.66 -0.75
is -0.28 -0.03 -0.33 0.67 1.00 0.80 -0.78
oi -0.24 -0.04 -0.35 0.66 0.80 1.00 -0.78
scs_ts 0.71 0.45 0.66 -0.75 -0.78 -0.78 1.00
In the COPE scale only two factors are identified.
all_items$pos_reinterpretation <- with(all_items, cope_1 + cope_29 + cope_38 + cope_59)
all_items$mental_disengagement <- with(all_items, cope_2 + cope_16 + cope_31 + cope_43)
all_items$venting <- with(all_items, cope_3 + cope_17 + cope_28 + cope_46)
all_items$seeking_instrumental_support <- with(all_items, cope_4 + cope_14 + cope_30 + cope_45)
all_items$active_coping <- with(all_items, cope_5 + cope_25 + cope_47 + cope_58)
all_items$denial <- with(all_items, cope_6 + cope_27 + cope_40 + cope_57)
all_items$religion <- with(all_items, cope_7 + cope_18 + cope_48 + cope_60)
all_items$humor <- with(all_items, cope_8 + cope_20 + cope_36 + cope_50)
all_items$behavioral_disengagement <- with(all_items, cope_9 + cope_24 + cope_37 + cope_51)
all_items$restraint <- with(all_items, cope_10 + cope_22 + cope_41 + cope_49)
all_items$seeking_emotional_support <- with(all_items, cope_11 + cope_23 + cope_34 + cope_52)
all_items$substance_use <- with(all_items, cope_12 + cope_26 + cope_35 + cope_53)
all_items$acceptance <- with(all_items, cope_13 + cope_21 + cope_44 + cope_54)
all_items$suppr_competing_activities <- with(all_items, cope_15 + cope_33 + cope_42 + cope_55)
all_items$planning <- with(all_items, cope_19 + cope_32 + cope_39 + cope_56)
Create COPE sub-scales scores using all items – note that SEM analyses suggest to drop some of the items.
all_items$active_coping <- with(
all_items, pos_reinterpretation + active_coping +
suppr_competing_activities + planning + restraint +
seeking_instrumental_support + acceptance
)
all_items$avoidance_coping <- with(
all_items, mental_disengagement + denial + humor +
behavioral_disengagement + substance_use + religion
)
all_items$soc_emo_coping <- with(
all_items, seeking_instrumental_support +
seeking_emotional_support + venting
)
plot(density(all_items$scs_ts))
fit_1 <- brm(
bf(
scs_ts ~ is_rescue_worker,
sigma ~ is_rescue_worker
),
family = student(),
backend = "cmdstanr",
data = all_items
)
pp_check(fit_1)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
me <- conditional_effects(
fit_1, "is_rescue_worker"
)
plot(me, points = FALSE)
summary(fit_1)
Family: student
Links: mu = identity; sigma = log; nu = identity
Formula: scs_ts ~ is_rescue_worker
sigma ~ is_rescue_worker
Data: all_items (Number of observations: 1068)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 74.11 1.02 72.15 76.10 1.00 3862 2674
sigma_Intercept 2.88 0.04 2.80 2.97 1.00 4140 2724
is_rescue_workerSi 7.49 1.19 5.17 9.78 1.00 4050 2963
sigma_is_rescue_workerSi -0.13 0.05 -0.23 -0.03 1.00 4027 2960
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
nu 38.40 15.75 16.97 77.94 1.00 4228 3115
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
BFt <- BayesFactor::ttestBF(
all_items$scs_ts[all_items$is_rescue_worker == "Si"],
all_items$scs_ts[all_items$is_rescue_worker == "No"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
------------------------
0.42 | [0.29, 0.56]
- Estimated using pooled SD.
rw_df <- all_items |>
dplyr::filter(is_rescue_worker == "Si")
rw_df <- rw_df %>%
mutate(job_qualification = case_when(
job_qualification == "non_rescue_worker" ~ "team_member",
TRUE ~ job_qualification))
lpa_scales <- c(
"is_rescue_worker",
"neuroticism", "extraversion", "openness", "agreeableness", "conscientiousness",
"active_coping", "avoidance_coping", "soc_emo_coping",
"iesr_ts",
# "avoiding", "intrusivity", "hyperarousal",
# "sk", "ch", "mi", "sj", "is", "oi",
# "pos_sc",
# "neg_sc",
# "ts_sc",
"mpss_tot"
# "ptgi_total_score"
# "relating_to_others",
# "new_possibilities",
# "personal_strength",
# "appreciation_of_life",
# "spirituality"
)
lpa_rw_df <- subset(rw_df, select=lpa_scales)
# lpa_rw_df <- subset(all_items, select=lpa_scales) |>
# dplyr::filter(is_rescue_worker == "Si")
lpa_rw_df$is_rescue_worker <- NULL
head(lpa_rw_df)
lpa_rw_df %>%
scale() %>%
estimate_profiles(1:6,
variances = c("equal", "varying"),
covariances = c("zero", "varying")
#package = "MplusAutomation"
)
tidyLPA analysis using mclust:
Model Classes AIC BIC Entropy prob_min prob_max n_min n_max BLRT_p
1 1 21342.45 21434.88 1.00 1.00 1.00 1.00 1.00
1 2 20740.50 20883.76 0.71 0.87 0.94 0.36 0.64 0.01
1 3 20459.15 20653.25 0.75 0.82 0.93 0.18 0.57 0.01
1 4 20314.46 20559.39 0.72 0.84 0.84 0.11 0.41 0.01
1 5 20232.20 20527.97 0.71 0.78 0.84 0.09 0.32 0.01
1 6 20161.53 20508.13 0.74 0.72 0.89 0.05 0.32 0.01
6 1 19882.54 20182.93 1.00 1.00 1.00 1.00 1.00
6 2 19509.14 20114.55 0.70 0.89 0.94 0.44 0.56 0.01
6 3 19415.46 20325.88 0.73 0.87 0.90 0.26 0.38 0.01
6 4 19396.37 20611.80 0.72 0.82 0.89 0.23 0.29 0.19
6 5 19359.41 20879.85 0.77 0.84 0.90 0.10 0.28 0.07
6 6 19321.00 21146.46 0.81 0.81 0.99 0.08 0.23 0.04
lpa_rw_df %>%
scale() %>%
estimate_profiles(1:6,
variances = c("equal", "varying"),
covariances = c("zero", "varying")
# package = "MplusAutomation"
) %>%
compare_solutions(statistics = c("AIC", "BIC"))
m2 <- lpa_rw_df %>%
scale() %>%
estimate_profiles(2,
variances = "varying",
covariances = "varying",
package = "MplusAutomation"
)
m2_plot <- lpa_rw_df %>%
scale() %>%
estimate_profiles(2,
variances = "varying",
covariances = "varying",
package = "MplusAutomation"
) %>%
plot_profiles(add_line = TRUE, rawdata= FALSE, bw = FALSE)
Profile 2: dysfunctional Profile 1: adaptive
get_estimates(m2)
out <- get_data(m2)
lpa_rw_df$lpa_class <- out$Class
table(
lpa_rw_df$lpa_class
)
1 2
409 342
table(
lpa_rw_df$lpa_class, rw_df$job_qualification
)
driver team_leader team_member
1 96 161 152
2 59 138 145
lpa_rw_df$class <- factor(lpa_rw_df$lpa_class)
summary(lpa_rw_df$class)
1 2
409 342
rw_df$class <- lpa_rw_df$class
rw_df$profile <- lpa_rw_df$class
rw_df$profile <- ifelse(
rw_df$profile == "2", "maladaptive", "adaptive"
)
rw_df$profile <- factor(rw_df$profile)
# Reorder the levels of the 'profile' factor
rw_df$profile <- factor(rw_df$profile, levels = c("maladaptive", "adaptive"))
m1 <- brm(
bf(scs_ts ~ profile),
family = skew_normal(),
data = rw_df,
init = 0.1,
backend = "cmdstanr",
adapt_delta = 0.9
)
pp_check(m1)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m1)
Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: scs_ts ~ profile
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 74.55 0.83 72.90 76.14 1.00 3443 2621
profileadaptive 12.64 1.15 10.39 14.92 1.00 3493 2699
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 14.90 0.38 14.19 15.68 1.00 4123 2947
alpha -0.70 0.59 -1.53 0.58 1.00 2187 3498
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m1, "profile"
)
plot(me, points = FALSE)
BFt <- BayesFactor::ttestBF(
rw_df$scs_ts[rw_df$class == 1],
rw_df$scs_ts[rw_df$class == 2],
paired = FALSE
)
effectsize(BFt, type = "d")
Cohen's d | 95% CI
------------------------
0.87 | [0.71, 1.02]
- Estimated using pooled SD.
m1 %>%
emmeans( ~ profile) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = profile, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
labs(x = "LPA Class", y = "SCS Score", title = "Rescue Workers") +
papaja::theme_apa() +
annotate("text", x = 1, y = 83, label = "Bayesian Cohen's d = 0.89\n 95% CI [0.73, 1.04]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
names(all_items)
m2 <- brm(
bf(sj ~ profile),
data = rw_df,
family = student,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m2)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
loo_m2 <- loo(m2)
plot(loo_m2)
summary(m2)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: sj ~ profile
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 17.03 0.23 16.58 17.48 1.00 3615 2904
profileadaptive -3.61 0.31 -4.20 -2.99 1.00 3754 2741
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 4.27 0.12 4.04 4.52 1.00 2939 2589
nu 46.35 16.93 21.54 86.78 1.00 3388 2924
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m2, "profile"
)
plot(me, points = FALSE)
emmeans(m2, specs = pairwise ~ profile)
$emmeans
profile emmean lower.HPD upper.HPD
maladaptive 17.0 16.6 17.5
adaptive 13.4 13.0 13.8
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
maladaptive - adaptive 3.61 3.01 4.21
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$sj[rw_df$profile == "adaptive"],
rw_df$sj[rw_df$profile == "maladaptive"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
--------------------------
-0.81 | [-0.97, -0.66]
- Estimated using pooled SD.
p2 <- m2 %>%
emmeans( ~ profile) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = profile, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "SCS Self-Judgment", title = "A") +
papaja::theme_apa() +
annotate("text", x = 1, y = 14.5, label = "Bayesian Cohen's d = 0.82\n 95% CI [0.67, 0.97]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p2
m3 <- brm(
bf(is ~ profile),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.95
)
pp_check(m3)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
loo_m3 <- loo(m3)
plot(loo_m3)
summary(m3)
Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: is ~ profile
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 12.29 0.25 11.78 12.80 1.00 2470 2490
profileadaptive -3.38 0.34 -4.05 -2.69 1.00 2358 2099
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 4.06 0.12 3.84 4.32 1.00 2494 2593
alpha 2.31 0.65 1.13 3.65 1.00 1915 2113
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m3, "profile"
)
plot(me, points = FALSE)
emmeans(m3, specs = pairwise ~ profile)
$emmeans
profile emmean lower.HPD upper.HPD
maladaptive 12.29 11.80 12.81
adaptive 8.91 8.49 9.28
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
maladaptive - adaptive 3.38 2.74 4.09
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$is[rw_df$profile == "adaptive"],
rw_df$is[rw_df$profile == "maladaptive"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
--------------------------
-0.93 | [-1.08, -0.78]
- Estimated using pooled SD.
p3 <- m3 %>%
emmeans( ~ profile) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = profile, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "SCS Isolation", title = "B") +
papaja::theme_apa() +
annotate("text", x = 1, y = 10.2, label = "Bayesian Cohen's d = 0.94\n 95% CI [0.79, 1.10]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p3
m4 <- brm(
bf(oi ~ profile),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m4)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
loo_m4 <- loo(m4)
plot(loo_m4)
summary(m4)
Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: oi ~ profile
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 10.96 0.24 10.48 11.42 1.00 1492 1549
profileadaptive -2.51 0.34 -3.16 -1.85 1.00 1366 1537
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 3.67 0.12 3.44 3.92 1.00 1369 1758
alpha 3.59 0.95 2.16 5.73 1.00 1253 1261
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m4, "profile"
)
plot(me, points = FALSE)
emmeans(m4, specs = pairwise ~ profile)
$emmeans
profile emmean lower.HPD upper.HPD
maladaptive 10.97 10.48 11.42
adaptive 8.44 8.04 8.79
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
maladaptive - adaptive 2.52 1.86 3.16
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$oi[rw_df$profile == "adaptive"],
rw_df$oi[rw_df$profile == "maladaptive"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
--------------------------
-0.96 | [-1.11, -0.81]
- Estimated using pooled SD.
p4 <- m4 %>%
emmeans( ~ profile) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = profile, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "SCS Over-Identification", title = "C") +
papaja::theme_apa() +
annotate("text", x = 1, y = 9.4, label = "Bayesian Cohen's d = 0.97\n 95% CI [0.82, 1.12]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p4
m5 <- brm(
bf(sk ~ profile),
family = student(),
data = rw_df,
backend = "cmdstanr"
)
pp_check(m5)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
loo_m5 <- loo(m5)
plot(loo_m5)
summary(m5)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: sk ~ profile
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 13.05 0.23 12.60 13.52 1.00 3991 3028
profileadaptive 1.25 0.31 0.63 1.86 1.00 4251 2936
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 4.21 0.12 3.98 4.45 1.00 3510 2614
nu 41.31 16.98 17.09 83.01 1.00 3245 2884
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m5, "profile"
)
plot(me, points = FALSE)
emmeans(m5, specs = pairwise ~ profile)
$emmeans
profile emmean lower.HPD upper.HPD
maladaptive 13.1 12.6 13.5
adaptive 14.3 13.9 14.7
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
maladaptive - adaptive -1.26 -1.85 -0.625
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$sk[rw_df$profile == "adaptive"],
rw_df$sk[rw_df$profile == "maladaptive"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
------------------------
0.28 | [0.13, 0.42]
- Estimated using pooled SD.
p5 <- m5 %>%
emmeans( ~ profile) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = profile, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "SCS Self-Kindness", title = "D") +
papaja::theme_apa() +
annotate("text", x = 1, y = 14.5, label = "Bayesian Cohen's d = 0.28\n 95% CI [0.14, 0.43]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p5
m6 <- brm(
bf(ch ~ profile),
family = student(),
data = rw_df,
backend = "cmdstanr"
)
pp_check(m6)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
loo_m6 <- loo(m6)
plot(loo_m6)
summary(m6)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: ch ~ profile
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 11.59 0.18 11.23 11.96 1.00 4291 3143
profileadaptive 0.01 0.24 -0.47 0.48 1.00 4425 3056
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 3.28 0.09 3.10 3.46 1.00 4625 2797
nu 43.03 16.67 19.04 83.45 1.00 4370 3134
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m6, "profile"
)
plot(me, points = FALSE)
BFt <- BayesFactor::ttestBF(
rw_df$ch[rw_df$profile == "adaptive"],
rw_df$ch[rw_df$profile == "maladaptive"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
-------------------------
-5.07e-04 | [-0.14, 0.14]
- Estimated using pooled SD.
p6 <- m6 %>%
emmeans( ~ profile) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = profile, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "SCS Common-Humanity", title = "E") +
papaja::theme_apa() +
annotate("text", x = 1.5, y = 12.2, label = "Bayesian Cohen's d = 0.00\n 95% CI [-0.14, 0.14]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p6
m7 <- brm(
bf(mi ~ profile),
family = student(),
data = rw_df,
backend = "cmdstanr"
)
Compiling Stan program...
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
Start sampling
Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: student_t_lpdf: Scale parameter is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 75, column 4 to column 48)
Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 1
Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: student_t_lpdf: Scale parameter is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 75, column 4 to column 48)
Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 1
Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 2 Exception: gamma_lpdf: Random variable is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 67, column 2 to column 66)
Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 2
Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 2 Exception: student_t_lpdf: Scale parameter is 0, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 75, column 4 to column 48)
Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 2
Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 3 Exception: gamma_lpdf: Random variable is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 67, column 2 to column 66)
Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 3
Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 3 Exception: student_t_lpdf: Scale parameter is 0, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 75, column 4 to column 48)
Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 3
Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 4 Exception: gamma_lpdf: Random variable is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 67, column 2 to column 66)
Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 4
pp_check(m7)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
loo_m7 <- loo(m7)
plot(loo_m7)
summary(m7)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: mi ~ profile
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 12.74 0.17 12.41 13.08 1.01 3843 2838
profileadaptive 1.01 0.23 0.58 1.46 1.00 3804 2789
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 3.03 0.09 2.86 3.21 1.00 3965 2848
nu 36.26 15.13 15.29 72.29 1.00 3682 3429
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
emmeans(m7, specs = pairwise ~ profile)
$emmeans
profile emmean lower.HPD upper.HPD
maladaptive 12.7 12.4 13.1
adaptive 13.8 13.5 14.1
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
maladaptive - adaptive -1.01 -1.44 -0.557
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$mi[rw_df$profile == "adaptive"],
rw_df$mi[rw_df$profile == "maladaptive"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
------------------------
0.32 | [0.18, 0.46]
- Estimated using pooled SD.
me <- conditional_effects(
m7, "profile"
)
plot(me, points = FALSE)
p7 <- m7 %>%
emmeans( ~ profile) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = profile, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "SCS Mindfulness", title = "F") +
papaja::theme_apa() +
annotate("text", x = 1, y = 13.8, label = "Bayesian Cohen's d = 0.32\n 95% CI [0.17, 0.46]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p7
fig_scs <- (p2 | p3 | p4) /
(p5 | p6 | p7)
out <- fig_scs + plot_annotation(
title = 'SCS Subscales as a Function of LPA Profile'
# subtitle = 'Rescue Workers group'
# caption = 'Disclaimer: None of these plots are insightful'
)
ggsave("scs_subscales_lpa.pdf", width = 35, height = 20, units = "cm")
print(out)
m10 <- brm(
bf(ies_ts ~ class),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m10)
summary(m10)
p10 <- m10 %>%
emmeans( ~ class) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = class, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "Impact of Event Scale - Revised (IES-R)") +
# papaja::theme_apa() +
annotate("text", x = 1, y = 17, label = "Bayesian Cohen's d = 1.34\n 95% CI [1.18, 1.50]")
p10
BFt <- BayesFactor::ttestBF(
rw_df$ies_ts[rw_df$class == 1],
rw_df$ies_ts[rw_df$class == 2],
paired = FALSE
)
effectsize(BFt)
emmeans(m10 , specs = pairwise ~ class)
m11 <- brm(
bf(ptgi_total_score | trunc(lb = 0) ~ class),
family = student(),
data = rw_df,
backend = "cmdstanr"
)
pp_check(m11)
summary(m11)
rw_df$job_qualification <- ifelse(
rw_df$job_qualification == "non_rescue_worker", "team_member",
rw_df$job_qualification
)
m12 <- brm(
bf(ies_ts ~ job_qualification),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m12)
summary(m12)
me <- conditional_effects(
m12, "job_qualification"
)
plot(me, points = FALSE)
emmeans(m12 , specs = pairwise ~ job_qualification)
m13 <- brm(
bf(scs_ts ~ job_qualification * class),
family = gaussian(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m13)
me <- conditional_effects(
m13, "job_qualification:class"
)
plot(me, points = FALSE)
summary(m13)
emmeans(m13 , specs = pairwise ~ job_qualification*class)
m14 <- brm(
bf(ptgi_total_score ~ ies_total_score | (1 | comme\)),
family = gaussian(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
summary(m14)
bayes_R2(m14)
conditional_effects(m14, "is:class")
mydf <- data.frame(
scs = scale(rw_df$scs_ts),
class = ifelse(rw_df$class == 1, 0.0, 1.0),
ptgi = scale(rw_df$ptgi_total_score),
psc = scale(rw_df$sk + rw_df$ch + rw_df$mi),
nsc = scale(rw_df$sj + rw_df$oi + rw_df$is),
commettee = rw_df$red_cross_commeetee_location,
id = 1:nrow(rw_df)
)
mydf <- mydf[complete.cases(mydf), ]
# Impute NAs with the mode.
# mydf$commettee[is.na(mydf$commettee)] <- "Comitato di Groane"
# summary(factor(mydf$commettee))
f1 <- bf(scs ~ class, family = skew_normal())
f2 <- bf(ptgi ~ scs + class, family = skew_normal())
mod <- brm(
f1 + f2 + set_rescor(FALSE),
data = mydf,
cores = 4,
refresh = 0,
init = 0.1,
backend = "cmdstanr",
adapt_delta = 0.99
)
bayestestR::mediation(mod, mediator = "scs", ci = 0.95, method = "SPI")
pp_check(mod, resp = "ptgi")
pp_check(mod, resp = "scs")
summary(mod)
fit <- mlm(d = mydf,
id = "commettee",
x = "class",
m = "scs",
y = "ptgi",
iter = 2000,
cores = 4)
mlm_path_plot(fit, level = .95, text = T,
xlab = "Resilience\nProfile",
mlab = "Self\nCompassion",
ylab = "PTG", digits = 2)
Treatment: class Mediator : scs Response : ptgi
Direct Effect (ADE) | -0.300 | [-0.464, -0.140] Indirect Effect (ACME) | 0.121 | [ 0.053, 0.187] Mediator Effect | 0.154 | [ 0.069, 0.230] Total Effect | -0.178 | [-0.325, -0.027]
m16 <- brm(
bf(ptgi_total_score ~ class),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m16)
summary(m16)
bayes_R2(m16)
m17 <- brm(
bf(ptgi_total_score ~ is_rescue_worker),
family = skew_normal(),
data = all_items,
backend = "cmdstanr",
adapt_delta = 0.99
)
summary(m17)
rw_df$rate_of_activity_num <- as.integer(rw_df$rate_of_activity_num)
temp <- rw_df[1:746, ]
temp$activity <- cut(
temp$rate_of_activity_num,
breaks = c(-1, 0, 1, 2),
labels = c("Low", "Medium", "High"),
include.lowest = TRUE,
ordered_result = TRUE
)
m19 <- brm(
bf(activity ~ class),
family=cumulative("logit"),
data = temp,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m19)
summary(m19)
m20 <- brm(
bf(class ~ last_training_num * scs_ts),
family=bernoulli(),
data = temp,
backend = "cmdstanr",
adapt_delta = 0.99
)
summary(m20)
conditional_effects(m20, "scs_ts:last_training_num")
table(rw_df$class, rw_df$last_training_num)